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A general method is presented for determining the dynamic 
torsional/axial response of linear structures composed of either 
tapered bars or shafts to transient excitations. The method consists 
of formulating and solving the dynamic problem in the Laplace 
transform domain by the finite element method and obtaining the 
response by a numerical inversion of the transformed solution. The 
derivation of the torsional and axial stiffness matrices is based on 
the exact solution of the transformed governing equation of motion, 
and it consequently leads to the exact solution of the problem. The 
solution permits treatment of the most practical cases of linear 
tapered bars and shafts, and employs modelling of structures with 
only one element per member which reduces the number of degrees 
of freedom involved. The effects of external viscous or internal 
viscoelastic damping are also taken into account. 


INTRODUCTION 

The static dynamic and stability analysis of nonuniform structures composed of 
tapered beams and/or bars has attracted considerable attention (Chu et. al. 
1970; Kounadis 1975; Sato 1980). A thorough presentation of developments 
pertinent to the dynamic behavior of tapered beams/bars has been presented 
by Kolousek (1973). Lately, GangaRao and Spyrakos (1986) determined the 
static and dynamic response of tapered flexural/axial members through an 
analytical technique applicable to the wide class of initial-boundary value 
problems governed by linear differential operators with variable coefficients. 
Besides analytical methods restricted to limited cases due to the involved 
equations of motion and the associated conditions, numerical methods such as 
the Finite Difference Method (FDM) (Liable 1985) and especially the Finite 
Element Method (FEM), have been successfully employed (Gallagher et al. 1970; 
Rough et al. 1979). The FEM appears to be more popular than the FDM since 
it presents several organizational advantages and handles boundary conditions 
easier. Use of the FEM has been primarily based on the approximate lumped 
or consistent mass representation and on displacement functions which are 
solutions of the static governing equations (Beaufait et al. 1970; Gupta 1985). 
Tapered members are considered as an assembly of uniform elements with 
known stiffnesses which are super-imposed to construct the stiffness of the 
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member. This stepped representation requires a relatively large number of 

elements to accurately determine the dynamic response. In the case of 
linearly tapered members, an alternative approach would be the use of exact 
stiffness matrices developed from the solution of the static gpverning equation 
of axial/flexural deformation (Just 1977; Holzer 1986). Recently, Banerjee and 
Williams (1986) developed exact dynamic stiffness matrices for the axial, 
torsional, and flexural vibration of tapered beams to harmonically varying 
forces. The approximate as well as the exact stiffness matrices developed by 
Banerjee can be used in a conventional modal analysis formulation to provide 
the response of tapered structures. Such an analysis, however, requires 
prior determination of the natural frequencies and nodal shapes that can be 
obtained by solving the free vibration problem (Bathe 1982). Alternative 
highly accurate and efficient FEM formulations, based on transformed dynamic 
stiffness matrices, have been successfully employed by Spyrakos and Beskos, 
(1982) and Tamma et. al. (1987) for the dynamic analysis of frameworks 
modelled with uniform elements and subjected to general transient forces. In 
their analysis, the transformed dynamic stiffness matrices were developed 
through application of either Fourier or Laplace transform with respect to time 
on the equation of motion of a beam element. The structural response in the 
time domain is obtained from the transformed stiffness equation and a 
numerical inversion. Therefore, such an approach retains the advantages of 
the direct stiffness method eliminating the need for prior solution of an 
eigenvalue problem. 

In this paper, the dynamic response of structures composed of tapered bars 
and shafts to transient axial and/or torsional forces is determined. The 
formulation considers the most practical cases of cross-sections and types of 
taper, and includes effects of both external viscous and internal viscoelastic 
damping. The analysis employs the FEM with dynamic stiffness matrices 
expressed in the Laplace transform domain. The derivation of the stiffness 
matrices is based on the exact solution of the axial or torsional tapered 
element governing equations expressed in terms of Bessel functions. Thus, 
modelling of the structure requires only one element per member which 
reduces the number of degrees of freedom involved and simplifies the 
modelling of the configuration. Furthermore, evaluation of the response from 
the stiffness equation leads to the "exact" solution of the problem. A 
numerical Laplace transform based on Durbin’s algorithm (Durbin 1974) is then 
used to determine the structural response in the time domain. Durbin’s 
algorithm was chosen since it allows an efficient and accurate direct and 
inverse numerical Laplace transform of general forcing functions (Beskos et al. 
1983). 


FORMULATION OF THE PROBLEM 

Consider the general tapered bar element a-b with a straight centroidal 
axis and directions of the principal axes being the same for all cross sections 
as shown in Figure 1. The variation of the cross-sectional area A(x) and 
polar second moment of area J(x) may be represented as 


A(x) = A a (1 + r J) m 
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where L is the length of the element and A a , J a denote the values of the 
cross-sectional area and polar second moment of area, respectively, at the 
cross-section a in Figure 1. Given the geometrical properties of the element 
at the end sections (a and b), the positive constants r and m can be 
evaluated from the expressions 



Figure 1. Geometry and sign convention of a general bar/shaft 
element 


Even though the developments presented in the following sections are 
valid for any value of m from equation (2), special emphasis will be placed on 
the practical cases of linear taper with m=l and m=2. The case of m=l 
corresponds to rectangular and I-sections, while m=2 pertains to circular as 
well as I-sections (Gupta 1985). 

Axial Vibration 


The equation of motion for a small amplitude, free axial vibration of a 
linear elastic tapered bar (m=l) is 
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where u = u(x,t) is the axial displacement of the bar and E, p are the modulus 

of elasticity and the mass density of the bar, respectively. Expressing, 
A(x) in terms of ( = 1 + r ^ equation (3) takes the form 
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The Laplace transform y (s) of a function y(t) is defined by 


y(s) = J y(t) e st dt (5) 

o 

where s is, in general, a complex number. Application of the Laplace 
transform with respect to time on equation (4), under the assumption of zero 
initial conditions, yields 


f 2 u 
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, _ + m ? u - s* f 2 u = 0, (6) 

where primes indicate differentiation with respect to the spatial variable f. 
The general solution of equation (6) can be obtained on the basis of the 
procedure indicated by Myers (1971). The resulting expression contains Bessel 
functions of the second kind with complex kernels which are not readily 
applicable for a concise development of an element stiffness matrix. Thus, 
after some algebraic manipulations and use of properties of Bessel functions 
(Abramovitz et al. 1965), one can arrive at the following concise form of the 
general solution: 
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where C and C are constants and k = . 

1 2 Z 

Adapting as positive directions of the nodal displacements and forces 
the ones shown in Figure 1, the evaluation of the axial stiffness matrix for 
the bar element a-b can be obtained by relating the axial displacements at the 
nodes a and b to the axial forces 
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through the displacement function u(s) given by equation (7) (Spyrakos et al. 

1982). An entry k. ., through Laplace transformed stiffness matrix, is 

1 J 

defined as the transformed force at the ith degree of freedom due to a unit 
transformed displacement at the jth degree of freedom while all the other 
transformed displacements are zero. Thus, with the sign convention of Figure 
1, the following element transformed nodal force-displacement relationship 

in terms of the dynamic stiffness influence coefficients coefficients can 
be obtained 


f* 


► 



p 

* 

F t (s) 

► - 

k 

1 1 

k 

1 2 

4 

u 

i 


F (s) 
2 


k 

2 1 

k 

2 2 


u 

2 


0 4 



i 


ft 

4 


where 


364 



( 9 ) 


-H {I k (b)K n (a) + I n (a) \(h)) 


k = k = - L —r~ 

12 21 a(l+r) k 


= -H ( 1+r ) m {I (b) K. (a) + I, (a) K (b) > 
n K K n 


with A, a, b and n given by 

EA 

n = - | ( 1+m) , H = -jp s(p/E)' /2 , a = ^ (/»/ E) l/2 

b = ( l+r)a, B = I k (a) K^b) - I R (b) K^a) (10) 

For the case m = 2, the stiffness influence coefficients can be expressed in 
terms of hyperbolic functions through the relationships (Abramovitz et al . 
1965) 
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Thus, after some computational effort, equations (9) take the form 
k = -H (cosh(ar) + a ’ sinh(ar)} 

k = K = H (1 + r) (12) 
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It should be noted that the dynamic stiffness influence coefficients Djj 

presented by Beskos and Narayanan (1983) for a uniform bar element can be 
easily deduced from equations (12) and (13) for a=b and r tending to zero. 

Torsional Vibration 

The equation of motion for free tortional vibration of a linear elastic 
tapered shaft with circular cross-sections is 
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where <|>(x,t) is the angular displacement, G is the shear modulus, and C, which 
is equal to one for circular cross-sections represents the torsional rigidity of 
the cross-section. When C is given appropriate values, equation (14) can also 
be utilized to approximate the torsional response of a number of other cross- 
sections. Substituting J(x) given by equation (2) and expressing x in terms 
of 4, equation (14) results 
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Application of Laplace transform on equation (15) leads to 
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Observing the similarity between the equations of motion (6) and (16) and 
following the procedure employed for the treatment of the axial vibration, one 
can obtain the tortional stiffness equation 
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where the dynamic stiffness influence coefficients k. . can be determined 

-*■ j 

from equation (9) by replacing the variables m, a, b and H with t, «, § and D, 
respectively, given by 
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The positive directions of the nodal torsions and angular displacements are 
depicted in Figure 1. 

For the case m = 2, the dynamic stiffness influence coefficients k. . 

-L J 

can be expressed in terms of hyperbolic functions with the aid of expressions 
(11). Thus, 
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It is_of interest to note that the torsional dynamic stiffness influence 
coefficient Djj, which are presented by Beskos et al. (1983) for a uniform 
element, can be deduced as a particular case of the kjj given by equation (19) 
for a-fi and r tending to zero. 


EFFECT OF DAMPING 

Both internal and external viscous damping can be accounted for by the 
transformed dynamic stiffness influence coefficients. For reasons of simpli- 
city, the material of the bar is assumed to be a Kelvin solid obeying the 
constitutive law (Flugge 1967) 

a = W ( £ + f d! )’ (21) 

where a is the stress, s is the strain, W represents either the modulus of 
elasticity E or the shear modulus G, and f is the damping coefficient. 
Equation (21) in the Laplane transform domain takes the form 

<r = W (1 + f s ) £ (22) 


which implies that internal viscous damping can be considered by replacing W 
with W(l+fs) in equations ( 6 ) and (15), respectively. 

When external viscous damping is present, the additional damping force, 
R, is introduced in the equations of motion. Denoting with c the coefficient of 
damping, R can be expressed as 


R = 



(23) 


Application of Laplace transform on equation (17) yields 
R = - csu 


(24) 
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Expressions (22) and (24) indicate that equations (8) and (17) can account 
for combined external viscous and internal viscoelastic damping by replacing 
the variables a and a with a and , respectively, where 
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In contrast to the conventional way of accounting for damping as a 
percentage of the critical damping either in a mode superposition analysis or 
in a direct integration procedure, the present formulation allows the 
assignment of different damping properties for each individual structural 
member* As a result, the dynamic behavior of linear structures can be 
efficiently simulated in a more rational way. 

FORMULATION OF THE PROBLEM 

Once the dynamic stiffness coefficients are defined, the dynamic problem 
of a bar/shaft can be formulated in the following static-like form in the 
Laplace transform domain 


{F ( s ) } = [k(s)J {u(s>>, (26) 

where (F(s)} and (u(s)} represent the Laplace transformed axial/torsional 
dynamic load and displacement vectors, respectively. After the transformed 

boundary conditions are applied, (u(s)} can be obtained from equation (26) 
by a matrix inversion of the dynamic stiffness matrix for a sequence of values 
of s. Then the response (u(t)} in the time domain can be determined by a 
numerical inversion of the Laplace transformed displacement vector* The 
response {u(t)} is the exact solution of the dynamic problem, since the 
dynamic stiffness matrices have been developed from the exact solutions of the 
transformed equations of motion* 

The numerical algorithm adopted herein to invert the transformed response 
has been developed by Durbin (1974). It combines both finite Fourier cosine 
and sine transforms and operates with complex values of s. Thus, it is more 
time consuming than other algorithms operating with real data. Nevertheless, 
Durbin’s algorithm has been chosen since it provides higher accuracy than 
real data algorithms, a feature which is crucial in dynamic problems involving 
excitations of a transient time variation. 

The above formulation is based on the assumption of zero initial 
conditions. However, consideration of non-zero initial conditions does not 
present any difficulty. In this case, the Laplace transform of equation (4) 
yields 
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where u(x,o) and u(x,o) are the initial axial displacement and velocity of the 
axial element, respectively. Thus, the initial conditions in the Laplace 
transform domain can be represented by a load distributed along the length of 
the element. The distributed load can be converted to a vector of equivalent 
nodal forces (fj(s)}, i = 1,2, through standard finite element procedures 

(Davies 1980). 


NUMERICAL EXAMPLES 

This section presents the solutions of numerical examples in order to 
illustrate the method and demonstrate its merits. The numerical computations 
where performed on a IBM 3081-D computer. 


EXAMPLE 1 

Consider the structural system in figure 2 which consists of one tapered 
and one uniform bar with rectangular cross-sections having a constant width 
b. The numerical data pertaining to this system is L = 10 in (25.4 cm), b = 
1.0 in (2.54 

V»>4 



Fig. 2 Geometry and loading of the structural system of example 1 

cm), hi = h 2 = 0.5 in (1.27 cm), h 3 = 2.5 in (6.35 cm), p - 0.002 lb-sec 2 /in 4 
(0.0214 kg/cm 3 ), E = 10 7 lb/in 2 (6.89xl0 5 N/mm 2 ), and f 30 = 10 6 lb (4.448xl0 6 
N). The values of the subscript i = 1,2,3 denote the element nodes as shown 
in figure 2. With the aid of the kjj and the dynamic^ stiffness influence 
coefficients for a uniform bar (Beskos et al, 1983), D^j, the equilibrium 
equations in the frequency domain can be written in the form 
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Equation (29) solved for u (s) yields 
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with w = 


Evaluation of the Bessel functions appearing in the kij stiffness coefficients 
involve complex kernels with s ranging from very small to very large 

arguments. Thus, accurate evaluation of the kij requires use of appropriate 

assymptotic expansions of the modified Bessel functions I}^(s) and k]^(s) 

(Watson 1966). 

The response u 3 (t) in the time domain is obtained by a numerical 

inversion using Durbin’s algorithm and is plotted in figure 3. The total CPU 
time, including the formulation of equation (28), was only 0.14 secs. In order 
to establish the accuracy of the method, the u 3 (t) is also determined by the 
NASTRAN computer code using a mesh of twenty equal elements for the 
tapered bar and ten for the uniform member. The total CPU time required by 
NASTRAN was 26.87 secs. The present method required considerably less CPU 
time than NASTRAN, since it modelled the structure with one element per 
member. As shown in figure 3 the results obtained by NASTRAN and the 
present method are almost identical. 
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Fig. 3 Axial response u 3 (t) vs time of example 1 


EXAMPLE 2 

Consider the structural system of figure 4 that is composed of one 
tapered and one uniform shaft with circular cross-sections and subjected to a 
concentrated step torque of magnitude T 30 = 10 6 lb-in (11.29xl0 6 N-cm). The 
geometry of the structure is described by the parameters L = 10 in (25.4 cm), 
R, = R 2 = 0.3989 in (1.013 cm), R 3 = 0.8921 in (2.266 cm), q = 0.002 lb-sec 2 /in 4 
(0.0214 kg/cm 3 ) and E = 10 7 lb/in a (6.89x10* N/mm J )._ The torque -f 3 (s) 
acting at node 3 causes the torsional deformation ? 3 (s) which can be 
evaluating from 


♦ 3 (s) = -T a (s) 
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Fig. 4 Geometry and loading of the structural system of 
example 2 


In equations (32) the subscript 1 pertains to node 1 of the uniform shaft 
element. Figure 5 shows the angular response * 3 (t) in the time domain 
obtained by numerical inversion of + 3 (s). The same figure also portrays 

results obtained by NASTRAN for a discretization of forty elements for the 
tapered shaft and ten for the uniform member. The total CPU time required 
by NASTRAN was 38.23 secs, while the present method required only 0.09 secs. 
Evaluation of the angular response * 3 (t) by the present method required less 
computational time than the evaluation of the axial response u 3 (t) in the first 
example. This can be primarily attributed to the functional form of the kjj 
and Kjj stiffness coefficients of equations (29) and (31), respectively. The 
former are expressed in terms of Bessel functions, while the latter consists of 
hyperbolic functions. It should be mentioned that results obtained by 
NASTRAN for a twenty element discretization of the tapered member did not 
provide sufficient accuracy. 



Time (seconds) 

Fig. 5 Angular response $ 3 (t) vs time of example 2 
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CONCLUSIONS 


In this work, the exact dynamic stiffness matrices of viscoelastic tapered 
bar and shaft elements are developed. These matrices can be incorporated in 
a finite element formulation to determine the response of structural systems to 
dynamic forces of a transient time variation. The formulation is performed in 
the Laplace transform domain resulting in a static-like relationship between 
the force and displacement vectors. The dynamic response is then obtained in 
the frequency domain numerically, and is subsequently evaluated in the time 
domain by a numerical inversion. Although the geometries and loading of the 
example problems presented are simple, the present method is general and 
applies to complicated situations. 

Within the realm of assumptions and limitations of linear theories, the k ij 
and kij dynamic stiffness coefficients lead to the exact solution of the 
problem, since they have been developed from the exact solutions of the 
transformed equations of motion. Thus, results obtained by the present 
method can be used to compare the accuracy of other numerical methods such 
as conventional finite element and finite differences methods. 

Use of the kjj and kjj coefficients accounts for the inertia and stiffness 
properties of the system members accurately, through a modelling that 
requires only one element per member. This is a significant advantage of the 
proposed method over conventional finite element methods employing a lumped 
or a consistent mass representation. Further, the method does not require 
the evaluation of nodal shapes or eigenvectors. Any inaccuracy of results can 
be primarily attributed to the accuracy of the numerical inversion algorithm 
used. 


The kjj and kij coefficients permit consideration of different levels of 
external or internal viscoelastic damping at each one of the axial/torsional 
members, the supports and the joints. Thus, one can control the response 
through a more rational estimation of the damping attributed to the individual 
structural members. 
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